Setup

set.seed(234876)

library(MCMCglmm)
## Loading required package: Matrix
## Loading required package: coda
## Loading required package: ape
library(tidyverse)
## Loading tidyverse: ggplot2
## Loading tidyverse: tibble
## Loading tidyverse: tidyr
## Loading tidyverse: readr
## Loading tidyverse: purrr
## Loading tidyverse: dplyr
## Conflicts with tidy packages ----------------------------------------------
## expand(): tidyr, Matrix
## filter(): dplyr, stats
## lag():    dplyr, stats
library(parallel)
library(knitr)

source("h2life_load_data.R")
## Parsed with column specification:
## cols(
##   setDate = col_date(format = ""),
##   flipDate = col_date(format = ""),
##   days = col_integer(),
##   age = col_integer(),
##   NewAge = col_integer(),
##   fID = col_character(),
##   id = col_character(),
##   sireid = col_character(),
##   damid = col_character(),
##   repl = col_character(),
##   treat = col_character(),
##   NstartF = col_integer(),
##   status = col_integer()
## )
genet_corr <- tibble(model = character(),
                     HS_STD = numeric(),
                     LY_STD = numeric(),
                     HS_LY = numeric(),
                     n_eff = numeric())

iter <- 2.5e5
burnin <- 2e4
thinning <- 500
chains <- 12

rerun <- TRUE

MANOVA analysis

Poster for fly meeting

if (rerun) {
  HS <- h2life %>% 
    filter(treat == "HS") %>% rename(Age_HS = NewAge) %>% 
    as.data.frame()
  LY <- h2life %>% 
    filter(treat == "LY") %>% rename(Age_LY = NewAge) %>% 
    as.data.frame()
  STD <- h2life %>%
    filter(treat == "STD") %>% rename(Age_STD = NewAge) %>% 
    as.data.frame()
  
  h2life_mrg <- full_join(full_join(HS, LY), STD)
  
  prior1 <- list(R = list(V = diag(3) * 1.002, nu = 1.002),
                 G = list(G1 = list(V = diag(3) * 1.002, nu = 0.002)))
  
  prior2 <- list(R = list(V = diag(3) / 3, nu = 1.002),
                 G = list(G1 = list(V = diag(3) / 3, nu = 0.002)))
  
  prior3 <- list(R = list(V = diag(3) * 1.002, nu = 1),
                 G = list(G1 = list(V = diag(3) * 0.5, nu = 0.002)))
  
  priors <- list(prior1, prior2, prior3)
  prior_names <- c("Tri: V = diag(3) * 1.002, nu = 0.02",
                   "Tri: V = diag(3) / 3, nu = 0.02",
                   "Tri: V = diag(3) * 0.5, nu = 0.02")
  
  for (ii in 1:length(priors)) {
    prior <- priors[[ii]]
    fm <- mclapply(1:chains, function(i) {
      MCMCglmm(cbind(Age_STD, Age_HS, Age_LY) ~ trait - 1,
               random = ~ us(trait):animal,
               rcov = ~ idh(trait):units,
               data = h2life_mrg,
               prior = prior,
               pedigree = pedigree,
               family = rep("gaussian", 3),
               nitt = iter,
               burnin = burnin,
               thin = thinning)
    }, mc.cores = chains)
    outfile <- paste0("tri_model_prior", ii, ".Rda")
    save(fm, file = outfile)
    
    re <- lapply(fm, function(m) m$VCV)
    re <- do.call(mcmc.list, re)
    re <- as.mcmc(as.matrix(re))
    n_eff <- effectiveSize(re)
    
    # STD vs. HS
    HS_STD <- re[ , "traitAge_HS:traitAge_STD.animal"] /
      sqrt(re[ , "traitAge_STD:traitAge_STD.animal"] *
             re[ , "traitAge_HS:traitAge_HS.animal"])
    
    # STD vs. LY
    LY_STD <- re[ , "traitAge_LY:traitAge_STD.animal"] /
      sqrt(re[ , "traitAge_STD:traitAge_STD.animal"] *
             re[ , "traitAge_LY:traitAge_LY.animal"])
    
    # LY vs. HS
    HS_LY <- re[ , "traitAge_HS:traitAge_LY.animal"] /
      sqrt(re[ , "traitAge_LY:traitAge_LY.animal"] *
             re[ , "traitAge_HS:traitAge_HS.animal"])
    
    genet_corr <- bind_rows(genet_corr,
                            tibble(model = prior_names[[ii]],
                                   HS_STD = median(HS_STD),
                                   LY_STD = median(LY_STD),
                                   HS_LY = median(HS_LY),
                                   n_eff = mean(n_eff)))
  }
  
  write_csv(genet_corr, path = "Genetic_Correlations.csv")
}
## Joining, by = c("sire", "dam", "animal", "treat")
## Joining, by = c("sire", "dam", "animal", "treat")

Analyze model

load("tri_model_prior1.Rda")

fe <- lapply(fm, function(m) m$Sol)
fe <- do.call(mcmc.list, fe)
plot(fe, ask = FALSE)

plot(fe[, 1, drop = FALSE], ask = FALSE)

plot(fe[, 2, drop = FALSE], ask = FALSE)

plot(fe[, 3, drop = FALSE], ask = FALSE)

# Extract random effects, convert to mcmc.list
re <- lapply(fm, function(m) m$VCV)
re <- do.call(mcmc.list, re)

plot(re[, 1, drop = FALSE], ask = FALSE)

plot(re[, 2, drop = FALSE], ask = FALSE)

plot(re[, 3, drop = FALSE], ask = FALSE)

autocorr.diag(re)
##           traitAge_STD:traitAge_STD.animal traitAge_HS:traitAge_STD.animal
## Lag 0                         1.0000000000                     1.000000000
## Lag 500                       0.0842692033                     0.076989902
## Lag 2500                     -0.0178249648                    -0.018361472
## Lag 5000                      0.0151468155                    -0.001760693
## Lag 25000                     0.0008825723                    -0.007553334
##           traitAge_LY:traitAge_STD.animal traitAge_STD:traitAge_HS.animal
## Lag 0                         1.000000000                     1.000000000
## Lag 500                       0.069088642                     0.076989902
## Lag 2500                     -0.024071358                    -0.018361472
## Lag 5000                      0.009317922                    -0.001760693
## Lag 25000                     0.003487370                    -0.007553334
##           traitAge_HS:traitAge_HS.animal traitAge_LY:traitAge_HS.animal
## Lag 0                        1.000000000                    1.000000000
## Lag 500                      0.083773312                    0.083790353
## Lag 2500                    -0.004399008                    0.004611076
## Lag 5000                    -0.022509111                   -0.008720982
## Lag 25000                   -0.016121417                   -0.007318574
##           traitAge_STD:traitAge_LY.animal traitAge_HS:traitAge_LY.animal
## Lag 0                         1.000000000                    1.000000000
## Lag 500                       0.069088642                    0.083790353
## Lag 2500                     -0.024071358                    0.004611076
## Lag 5000                      0.009317922                   -0.008720982
## Lag 25000                     0.003487370                   -0.007318574
##           traitAge_LY:traitAge_LY.animal traitAge_STD.units
## Lag 0                        1.000000000        1.000000000
## Lag 500                      0.083268922        0.075959690
## Lag 2500                     0.005213926       -0.012102331
## Lag 5000                     0.007473910        0.015920811
## Lag 25000                   -0.014482780        0.004354346
##           traitAge_HS.units traitAge_LY.units
## Lag 0           1.000000000      1.0000000000
## Lag 500         0.075349367      0.0642032980
## Lag 2500        0.007898546     -0.0153543256
## Lag 5000       -0.024638971     -0.0001941856
## Lag 25000      -0.019878975     -0.0175330898
effectiveSize(re)
## traitAge_STD:traitAge_STD.animal  traitAge_HS:traitAge_STD.animal 
##                         5004.529                         5119.084 
##  traitAge_LY:traitAge_STD.animal  traitAge_STD:traitAge_HS.animal 
##                         5161.457                         5119.084 
##   traitAge_HS:traitAge_HS.animal   traitAge_LY:traitAge_HS.animal 
##                         4858.866                         4784.512 
##  traitAge_STD:traitAge_LY.animal   traitAge_HS:traitAge_LY.animal 
##                         5161.457                         4784.512 
##   traitAge_LY:traitAge_LY.animal               traitAge_STD.units 
##                         4899.712                         4807.313 
##                traitAge_HS.units                traitAge_LY.units 
##                         4849.351                         5019.077
# Concatenate samples
re <- as.mcmc(as.matrix(re))

head(re)
## Markov Chain Monte Carlo (MCMC) output:
## Start = 1 
## End = 7 
## Thinning interval = 1 
##      traitAge_STD:traitAge_STD.animal traitAge_HS:traitAge_STD.animal
## [1,]                        0.4146686                      0.12177726
## [2,]                        0.4518303                      0.11549431
## [3,]                        0.3258092                      0.09183121
## [4,]                        0.3272157                      0.05843721
## [5,]                        0.6759533                      0.26814659
## [6,]                        0.5333801                      0.23382353
## [7,]                        0.3621810                      0.15351107
##      traitAge_LY:traitAge_STD.animal traitAge_STD:traitAge_HS.animal
## [1,]                      0.10732055                      0.12177726
## [2,]                      0.05899287                      0.11549431
## [3,]                      0.18131683                      0.09183121
## [4,]                      0.15226261                      0.05843721
## [5,]                      0.25923669                      0.26814659
## [6,]                      0.24500887                      0.23382353
## [7,]                      0.15202631                      0.15351107
##      traitAge_HS:traitAge_HS.animal traitAge_LY:traitAge_HS.animal
## [1,]                      0.3561362                     0.10024041
## [2,]                      0.2730102                     0.05165598
## [3,]                      0.3557770                     0.13505231
## [4,]                      0.1921802                     0.09424014
## [5,]                      0.3621698                     0.10099553
## [6,]                      0.4371556                     0.30262471
## [7,]                      0.4382596                     0.24309716
##      traitAge_STD:traitAge_LY.animal traitAge_HS:traitAge_LY.animal
## [1,]                      0.10732055                     0.10024041
## [2,]                      0.05899287                     0.05165598
## [3,]                      0.18131683                     0.13505231
## [4,]                      0.15226261                     0.09424014
## [5,]                      0.25923669                     0.10099553
## [6,]                      0.24500887                     0.30262471
## [7,]                      0.15202631                     0.24309716
##      traitAge_LY:traitAge_LY.animal traitAge_STD.units traitAge_HS.units
## [1,]                      0.1728861          0.5121829         0.5017728
## [2,]                      0.2074450          0.5128448         0.5476254
## [3,]                      0.3866927          0.5540212         0.5016486
## [4,]                      0.2524715          0.5582436         0.5875449
## [5,]                      0.2146196          0.3887915         0.4857655
## [6,]                      0.3120685          0.4636957         0.4390513
## [7,]                      0.2938009          0.5428062         0.4631711
##      traitAge_LY.units
## [1,]         0.6115914
## [2,]         0.5809705
## [3,]         0.4701095
## [4,]         0.5232854
## [5,]         0.5458642
## [6,]         0.5432736
## [7,]         0.5280122
# STD vs. HS
plot(re[ , "traitAge_HS:traitAge_STD.animal"])

plot(re[ , "traitAge_STD:traitAge_STD.animal"])

HS_STD <- re[ , "traitAge_HS:traitAge_STD.animal"] /
  sqrt(re[ , "traitAge_STD:traitAge_STD.animal"] *
         re[ , "traitAge_HS:traitAge_HS.animal"])
plot(HS_STD)

median(HS_STD)
## [1] 0.2917385
HPDinterval(HS_STD)
##            lower     upper
## var1 0.006316709 0.5387122
## attr(,"Probability")
## [1] 0.95
# STD vs. LY
LY_STD <- re[ , "traitAge_LY:traitAge_STD.animal"] /
  sqrt(re[ , "traitAge_STD:traitAge_STD.animal"] *
         re[ , "traitAge_LY:traitAge_LY.animal"])
plot(LY_STD)

median(LY_STD)
## [1] 0.4928118
HPDinterval(LY_STD)
##          lower     upper
## var1 0.2525297 0.7072765
## attr(,"Probability")
## [1] 0.95
# LY vs. HS
HS_LY <- re[ , "traitAge_HS:traitAge_LY.animal"] /
  sqrt(re[ , "traitAge_LY:traitAge_LY.animal"] *
         re[ , "traitAge_HS:traitAge_HS.animal"])
plot(HS_LY)

median(HS_LY)
## [1] 0.348081
HPDinterval(HS_LY)
##           lower     upper
## var1 0.06371274 0.5937976
## attr(,"Probability")
## [1] 0.95

Plot for poster

library(tidyverse)
library(cowplot)

M <- data_frame(`HS vs. STD` = as.numeric(HS_STD),
                `LY vs. STD` = as.numeric(LY_STD),
                `HS vs. LY` = as.numeric(HS_LY))

colMeans(M)
## HS vs. STD LY vs. STD  HS vs. LY 
##  0.2864158  0.4844311  0.3401840
M %>% gather(Comparison, value) %>%
  ggplot(aes(value, color = Comparison)) +
  geom_line(stat = "density", size = 2) +
  labs(x = "Genetic Correlation", y = "Density") +
  theme(legend.position = c(0.15, 0.85),
        text = element_text(size = 24),
        legend.text = element_text(size = 18))

ggsave(last_plot(), file = "Genetic_Correlation_Plot.pdf",
       width = 8, height = 6)

Model following Ingleby et al. 2013

Setup data as above

V_nu_to_a_b <- function(V, nu) {
  require(tidyverse)
  require(invgamma)
  
  a <- nu / 2
  b <- (nu * V) / 2
  
  p <- tibble(
    x = seq(0, 10, length = 200),
    y = dinvgamma(x, a, b)) %>% 
    ggplot(aes(x, y)) +
    geom_line()
  print(p)
  
  return(list(alpha = a, beta = b))
}
V_nu_to_a_b(1, 0.002)
## Loading required package: invgamma
## Warning: Removed 1 rows containing missing values (geom_path).

## $alpha
## [1] 0.001
## 
## $beta
## [1] 0.001
genet_corr <- read_csv("Genetic_Correlations.csv")
## Parsed with column specification:
## cols(
##   model = col_character(),
##   HS_STD = col_double(),
##   LY_STD = col_double(),
##   HS_LY = col_double(),
##   n_eff = col_double()
## )
if (rerun) {
  
  prior1 <- list(R = list(V = 1, nu = 0.002),
                 G = list(G1 = list(V = diag(3) / 3, nu = 0.002)))
  
  prior2 <- list(R = list(V = 1, nu = 0.002),
                 G = list(G1 = list(V = diag(3) * 0.5 * 0.7, nu = 0.002)))
  
  M <- matrix(0.5 * 0.7, 3, 3)
  diag(M) <- 0.5
  prior3 <- list(R = list(V = 1, nu = 0.002),
                 G = list(G1 = list(V = M, nu = 0.002)))
  
  priors <- list(prior1, prior2, prior3)
  prior_names <- c("Sire: V = diag(3) / 3, nu = 0.002",
                   "Sire: V = diag(3) * 0.5 * 0.7, nu = 0.002",
                   "Sire: V = 0.5 diag, 0.35 off-diag, nu = 0.002")
  
  for (ii in 1:length(priors)) {
    prior <- priors[[ii]]
    
    fm <- mclapply(1:chains, function(i) {
      MCMCglmm(NewAge ~ treat,
               random = ~ us(treat):sire,
               data = h2life,
               prior = prior,
               family = "gaussian",
               nitt = iter,
               burnin = burnin,
               thin = thinning,
               verbose = TRUE)
    }, mc.cores = chains)
    
    outfile <- paste0("sire_model_prior", ii, ".Rda")
    save(fm, file = outfile)
    
    re <- lapply(fm, function(m) m$VCV)
    re <- do.call(mcmc.list, re)
    re <- as.mcmc(as.matrix(re))
    n_eff <- effectiveSize(re)
    
    # STD vs. HS
    HS_STD <- re[ , "treatHS:treatSTD.sire"] /
      sqrt(re[ , "treatHS:treatHS.sire"] *
             re[ , "treatSTD:treatSTD.sire"])
    median(HS_STD)
    
    # STD vs. LY
    LY_STD <- re[ , "treatSTD:treatLY.sire"] /
      sqrt(re[ , "treatSTD:treatSTD.sire"] *
             re[ , "treatLY:treatLY.sire"])
    median(LY_STD)
    
    # LY vs. HS
    HS_LY <- re[ , "treatHS:treatLY.sire"] /
      sqrt(re[ , "treatLY:treatLY.sire"] *
             re[ , "treatHS:treatHS.sire"])
    median(HS_LY)
    
    genet_corr <- bind_rows(genet_corr,
                            tibble(model = prior_names[[ii]],
                                   HS_STD = median(HS_STD),
                                   LY_STD = median(LY_STD),
                                   HS_LY = median(HS_LY),
                                   n_eff = mean(n_eff)))
    
  }
  write_csv(genet_corr, path = "Genetic_Correlations.csv")
  
  kable(genet_corr, digits = 3)
}
model HS_STD LY_STD HS_LY n_eff
Tri: V = diag(3) * 1.002, nu = 0.02 0.298 0.494 0.347 4776.086
Tri: V = diag(3) / 3, nu = 0.02 0.295 0.499 0.352 4820.542
Tri: V = diag(3) * 0.5, nu = 0.02 0.294 0.491 0.351 4812.753
Sire: V = diag(3) / 3, nu = 0.002 0.515 0.486 0.562 5436.460
Sire: V = diag(3) * 0.5 * 0.7, nu = 0.002 0.514 0.493 0.566 5499.585
Sire: V = 0.5 diag, 0.35 off-diag, nu = 0.002 0.520 0.489 0.566 5520.000

Analyze model

load("sire_model_prior1.Rda")

fe <- lapply(fm, function(m) m$Sol)
fe <- do.call(mcmc.list, fe)
plot(fe, ask = FALSE)

# Extract random effects, convert to mcmc.list
re <- lapply(fm, function(m) m$VCV)
re <- do.call(mcmc.list, re)

plot(re[, 1:3], ask = FALSE)

plot(re[, 4:6], ask = FALSE)

plot(re[, 7:9], ask = FALSE)

plot(re[, 10], ask = FALSE)

effectiveSize(re)
##   treatHS:treatHS.sire   treatLY:treatHS.sire  treatSTD:treatHS.sire 
##               5308.632               5520.000               5565.375 
##   treatHS:treatLY.sire   treatLY:treatLY.sire  treatSTD:treatLY.sire 
##               5520.000               5274.867               5358.924 
##  treatHS:treatSTD.sire  treatLY:treatSTD.sire treatSTD:treatSTD.sire 
##               5565.375               5358.924               5534.616 
##                  units 
##               5358.443
# Concatenate samples
re <- as.mcmc(as.matrix(re))
colnames(re)
##  [1] "treatHS:treatHS.sire"   "treatLY:treatHS.sire"  
##  [3] "treatSTD:treatHS.sire"  "treatHS:treatLY.sire"  
##  [5] "treatLY:treatLY.sire"   "treatSTD:treatLY.sire" 
##  [7] "treatHS:treatSTD.sire"  "treatLY:treatSTD.sire" 
##  [9] "treatSTD:treatSTD.sire" "units"
# ??? Heritability ????
HS <- re[, "treatHS:treatHS.sire"] / (re[, "treatHS:treatHS.sire"] + re[, "units"])
median(HS)
## [1] 0.1500825
LY <- re[, "treatLY:treatLY.sire"] / (re[, "treatLY:treatLY.sire"] + re[, "units"])
median(LY)
## [1] 0.09315043
STD <- re[, "treatSTD:treatSTD.sire"] / (re[, "treatSTD:treatSTD.sire"] + re[, "units"])
median(STD)
## [1] 0.1680397
# STD vs. HS
HS_STD <- re[ , "treatHS:treatSTD.sire"] /
  sqrt(re[ , "treatHS:treatHS.sire"] *
         re[ , "treatSTD:treatSTD.sire"])
plot(HS_STD)

median(HS_STD)
## [1] 0.5150163
HPDinterval(HS_STD)
##          lower     upper
## var1 0.1742536 0.8009837
## attr(,"Probability")
## [1] 0.95
# STD vs. LY
LY_STD <- re[ , "treatSTD:treatLY.sire"] /
  sqrt(re[ , "treatSTD:treatSTD.sire"] *
         re[ , "treatLY:treatLY.sire"])
plot(LY_STD)

median(LY_STD)
## [1] 0.486466
HPDinterval(LY_STD)
##          lower    upper
## var1 0.1232235 0.784221
## attr(,"Probability")
## [1] 0.95
# LY vs. HS
HS_LY <- re[ , "treatHS:treatLY.sire"] /
  sqrt(re[ , "treatLY:treatLY.sire"] *
         re[ , "treatHS:treatHS.sire"])
plot(HS_LY)

median(HS_LY)
## [1] 0.5620322
HPDinterval(HS_LY)
##          lower     upper
## var1 0.2035833 0.8322286
## attr(,"Probability")
## [1] 0.95